In seismic imaging, as well as in many other areas of seismic data processing, one often needs to compute integrals of the general form:
                              ∫          D                ⁢                  A          ⁢                                          ⁢                      ⅆ            V                                              (        1        )            
In equation (1), A is a function to be integrated, and dV is a volume element on D. The integration domain D is a k-dimensional sub-domain of n-dimensional space (k≦n).
The integration can be over, for example, source and receiver positions, midpoints only, midpoint and offsets or scattering angles etc. The integration domain typically is a two-dimensional domain, but it could be of higher dimension. The integrand A contains the seismic data, or part of the seismic data, and may also contain some theoretical terms. Integrals of this type occur in seismic data processing in, for example, migration algorithms (pre stack and post stack migration; time and depth migration; Kirchhoff, beam and wave equation migration), attenuation of multiple reflections in marine data, AVO analysis etc. Where k<n the topography of the acquisition geometry is taken into account—this happens when not only the x and y coordinates of the sources and receivers, but also their z coordinates, are known.
One particular problem in processing seismic data is that the seismic data, and hence the value of the function A, are generally known only at np n-dimensional points in the integration domain D. These points at which the seismic data are known will hereinafter be referred to as “data points”. Integrals having the form of equation (1) have traditionally been computed using a ‘binning’ algorithm as described in, for example, “Seismic data processing” by O. Yilmaz in Volume II, Society of Exploration Geophysicists, Tulsa, Okla. p 1019 (2001). A binning algorithm can be implemented in various ways. For example, one can evaluate equation (1) by using a summation in which all the data points are given equal weights:
                                          ∫            D                                                          ⁢                      A            ⁢                          ⅆ              V                                      ⁢                                  =                              Vol            ⁡                          (              D              )                                ⁢                                    ∑                              i                =                1                                            n                p                                      ⁢                                                  ⁢                                          A                i                            .                                                          (        2        )            
Vol(D) can be estimated roughly, but is often left out altogether in which case equation (2) reduces to:
                                          ∫            D                                                          ⁢                      A            ⁢                          ⅆ              V                                      ⁢                                  =                              ∑                          i              =              1                                      n              p                                ⁢                                          ⁢                      A            i                                              (        3        )            
Alternatively one can divide the integration domain D up in a large number of ‘bins’ (typically squares or rectangles in the case of a two-dimensional domain), after which the average value of the function A in each bin is computed (see e.g. Yilmaz, 2001). The integral is then evaluated as the sum of all these average values:
                                          ∫            D                                                          ⁢                      A            ⁢                          ⅆ              V                                      ⁢                                  =                              ∑            i                    ⁢                                          ⁢                                                    A                ~                            i                        ⁢                                          Vol                ⁡                                  (                                      b                    i                                    )                                            .                                                          (        4        )            
In equation (4) Ãi denotes the average of all measured values of A in the ith bin, Vol(bi) is the volume of the ith bin and the summation is over all the bins.
These prior art binning techniques can work reasonably well if the data points are distributed on a regular grid. However, the results are dependent on the binning grid size so that, even if the data points are distributed on a regular grid, errors can result if an unsatisfactory binning grid is applied. Furthermore, the data points acquired in a seismic survey are often not distributed over a regular grid—firstly, it is hard in practice to position the seismic sources and receivers exactly on a regular grid and, secondly, the actual optimal survey grid may be irregular rather than regular. If the data points are irregularly distributed over the integration domain then the prior art binning techniques are subject to further errors and become unreliable.